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Abstract 



Several relaxation approximations to partial differential equations have 
been recently proposed. Examples include conservation laws, ffamilton- 
Jacobi equations, convection-diffusion problems, gas dynamics problems. 
The present paper focuses onto diffusive relaxation schemes for the numer- 
^vj . ical approximation of nonlinear parabolic equations. These schemes are 

^ ' based on a suitable semilinear hyperbolic system with relaxation terms. 

^SJ , High order methods are obtained by coupling ENO and WENO schemes 

t~^ . for space discretization with IMEX schemes for time integration. Error es- 

lO ' timates and convergence analysis are developed for semidiscrete schemes 

\l , with numerical analysis for fully discrete relaxed schemes. Various nu- 

^■^ ' merical results in one and two dimensions illustrate the high accuracy and 

good properties of the proposed numerical schemes, also in the degenerate 
case. These schemes can be easily implemented on parallel computers and 
f^ ' applied to more general systems of nonlinear parabolic equations in two- 

"ti ' and three-dimensional cases. 

s 

j>I ; 1 Introduction 

^\ ' Relaxation approximations to nonlinear partial differential equations have been 

S . introduced ^ |21 I21L \'2'6\ 1251 1271 |2H1 on the basis of the replacement of the equa- 

tions with a suitable semilinear hyperbolic system with stiff relaxation terms. 
The idea can be explained by considering the case of a scalar conservation law 

du df{u) ^^ 
dt dx 



^Dipartimento di Matematica, Universita di Milano, Via Saldini 50, 1-20100 Milano, Italy. 

■f Dipartimento di Matematica, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 
Torino, Italy, gabriella.puppo@polito.it 

§This work was partially supported by the MIUR/PRIN2005 project "ModcUistica numer- 
ica per il calcolo scientifico ed applicazioni avanzate" . 



for which Jin and Xin proposed in |23| the following relaxation approximation 

du dv 
dt dx 

where e (relaxation time) is a small positive parameter and a is a positive 
constant satisfying 

for all u. For small values of e the previous system gives the following first order 
approximation of the original conservation law 

The generalization of the above model to multidimensional systems of conser- 
vation laws can be done in a natural way by adding more rate equations. One 
would expect that appropriate numerical schemes for the relaxation system yield 
accurate numerical approximations to the original equation or system when the 
relaxation rate e is sufficiently small. Numerically, the main advantage of solv- 
ing the relaxation model over the original conservation law lies in the simple 
linear structure of characteristic fields and the localized lower order term. In 
particular, the semilinear nature of the relaxation system gives a new way to 
develop numerical schemes that are simple, general and Riemann solver free 
18 20 23j . Several other relaxation approximation have been introduced re- 
cently. For example we mentioned here the work of Coquel and Perthame [TU] 
for real gas computation, the relaxation schemes of Jin et al. |19| for curvature- 
dependent front propagation, the relaxation approximation in the rapid granular 
flow |22j , the relaxation approximation and relaxation schemes for diffusion and 
convection-diffusion problems |21l 1251 l?7l 1281 . Moreover, there are strong and 
interesting links between the relaxation approximation and the kinetic approach 
to nonlinear transport equations, based upon analogies with the passage from 
the Boltzmann equation to fluid mechanics (see for example mUE]). 

The aim of this work is to analyze from both a theoretical and computational 
point of view relaxation schemes that approximate the following nonlinear de- 
generate parabolic problem 

^ = D/^(p{u)), rr e R^ i > 0, (1) 

with initial data m(x,0) = uq{x) E L^(R'*) and _D > is a diffusivity coefficient. 
As usual, we will assume p : M ^ M to be non-decreasing and Lipschitz continu- 
ous [321 ■ The equation is degenerate if p(0) = 0. We will also consider the same 
parabolic problem in a bounded domain O C W^ with Neumann boundary con- 
ditions. In the case p{u) — m™, to > 1 the previous equation is the porous media 
equation which describes the flow of a gas through a porous interface according 



to some constitutive relation like Darcy's law in order to link the velocity of the 
gas and its pressure. In this case the diffusion coefficient mu™'~^ vanishes at 
the points where u = and the governing parabolic equation degenerates there. 
The set of such points is called interface. Moreover the porous media equation 
can exhibit a finite speed of propagation for compactly supported initial data 
P). The influence of the degenerate diffusion terms make the dynamics of the 
interfaces difhcult to study from both the theoretical and the numerical point 
of view. Another interesting case corresponds to < m < 1 and it is referred to 
as the fast diffusion equation which appears, for example, in curvature-driven 
evolution and avalanches in sandpiles. In general the numerical analysis of equa- 
tion (Q) is difficult for at least two reasons: the appearance of singularities for 
compactly supported solutions and the growth of the size of the support as time 
increases [retention property). 

From the numerical viewpoint, an usual technique to approximate Q in- 
volves implicit discretization in time: it requires, at each time step, the dis- 
cretization of a nonlinear elliptic problem. However, when dealing with nonlin- 
ear problems one generally tries to linearize them in order to take advantage 
of efficient linear solvers. Linear approximation schemes based on the so-called 
nonlinear Chernoff's formula with a suitable relaxation parameter have been 
studied for example in jHl I29L I3UI ESj where also some energy error estimates 
have been investigated. Other linear approximation schemes have been intro- 
duced by Jager, Kacur and Handlovicova |17[I24| . More recently, different ap- 
proaches based on kinetic schemes for degenerate parabolic systems have been 
considered and analyzed by Aregba-DrioUet, Natalini and Tang in |5]. Other 
approaches were investigated in the work of Karlsen et al. ^H] based on a suit- 
able splitting technique with applications to more general hyperbolic-parabolic 
convection-diffusion equations. Finally, a new scheme based on the maximum 
principle and on a perturbation and regularization approach was proposed by 
Pop and Yong in |32,. 

Our approach is inspired by relaxation schemes where the nonlincarity inside 
the equation is replaced by a semilinearity. This reduction is carried out in 
order to obtain numerical schemes that are easy to implement, also for parallel 
computing, even in the multidimensional case and for more general and complex 
problems, like oil recovery problems |I2| . Moreover with our approach it is 
possible to improve numerical schemes by using high order methods and, in 
principle, different numerical approaches (finite volume, finite differences,. . . ). 
In particular, in this paper, in order to obtain high order methods, we couple 
ENO and WENO schemes for space discretization and IMEX schemes for time 
advancement. High order schemes may not reach their order of convergence due 
to the loss of regularity of the solution during the evolution. However they are 
interesting nevertheless for error reduction when the number of grid points is 
fixed or until discontinuities develop (both cases arise for example in nonlinear 
filtering in image analysis [SBJ). 

The paper is organized as follows. Section [5] is devoted to the introduction 
of our relaxation schemes. The stability and error estimates of the semidiscrete 
scheme is provided in Section 13 In Section 0] we consider the fully discrete 



relaxed scheme with nonhnear stabihty analysis, the numerical treatment of 
boundary conditions and the extension to the multi-dimensional case. Finally, 
the implementation of the method as well as the results of several numerical 
experiments are discussed in section |S1 

2 Relaxation approximation of nonlinear diffu- 
sion 



The main purpose of this work is to approximate solutions of a nonlinear, possi- 
bly degenerate, parabolic equation of the form ^, This framework is so general 
that it includes the porous medium equation p{u) ~ m'" with to > 1 ([3S| and 
references therein), non linear image processing [201, as well as a wide class of 
mildly nonlinear parabolic equations jl4| . 

The schemes proposed in the present work are based on the same idea at the 
basis of the well-known relaxation schemes for hyperbolic conservation laws \2'6\ . 
In the case of the nonlinear diffusion operator, an additional variable v{x, t) G W^ 
and a positive parameter e are introduced and the following relaxation system 
is obtained: 



du 
'dt 



— + div(w) = 



dv 
'dt 



— Vp(m) 



1^ 

'—V 

£ 



(1) 



Formally, in the small relaxation limit, £ ^ 0^, system (^ approximates to 
leading order equation JQ). In order to have bounded characteristic velocities 
and to avoid a singular differential operator as £ ^ 0+, a suitable parameter (p 
is introduced and JQ) can be rewritten as: 

+ div(w) = 

(2) 



dt 



dv 



1 



— + (^ Vp(m) = — V 

Ot £ 



V' 



- J ] Vp(u) 



Finally we remove the non linear term from the convective part, as in standard 
relaxation schemes, introducing a variable w{x, i) G M and rewriting the system 
as: 

— + div(z;) = 



Formally, as e 
recovered. 



dt 
dv 

dw 
'dt 
0+, w 



■ (p Ww = - 

- div('i/) = - 
-^ P{u), V 



1 



1 



D 



Vw 



(3) 



iw-p{u)) 



—Vp(u) and the original equation is 



In the previous system the parameter e has physical dimensions of time 
and represents the so-called relaxation time. Furthermore, w has the same 
dimensions as u, while each component of v has the dimension of u times a 
velocity; finally </? is a velocity. The inverse of e gives the rate at which v decays 
onto — Vp(w) in the evolution of the variable v governed by the stiff second 
equation of ^. 

Equations ^ form a semilinear hyperbolic system with a stiff source term. 
The characteristic velocities of the hyperbolic part are given by 0, ±(/j. The 
parameter ip allows to move the stiff terms —yp{u) to the right hand side, 
without losing the hyperbolicity of the system. 

We point out that degenerate parabolic equations often model physical sit- 
uations with free boundaries or discontinuities: we expect that schemes for 
hyperbolic systems will be able to reproduce faithfully these details of the so- 
lution. One of the main properties of Q consists in the semilinearity of the 
system, that is all the nonlinearities are in the (stiff) source terms, while the 
differential operator is linear. Hence, the solution of the convective part requires 
neither Riemann solvers nor the computation of the characteristic structure at 
each time step, since the eigenstructure of the system is constant in time. More- 
over, the relaxation approximation does not exploit the form of the nonlinear 
function p and hence it gives rise to a numerical scheme that, to a large extent, 
is independent of it, resulting in a very versatile tool. 

We also anticipate here that, in the relaxed case (i.e. e = 0), the stiff source 
terms can be integrated solving a system that is already in triangular form and 
then it does not require iterative solvers. 

3 The semidiscrete scheme 

System ^ can be written in the form: 

zt+div/(z) = ig(z), (1) 

where 

f{z)= $^w 9{z)=\ -v + {ip^e-D)Vw (2) 

p{u) — w J 

and $^ is the d x d identity matrix times the scalar ifP' . We start discretizing 
the system in time using, for simplicity, a uniform time step At. Let z^^(x) = 
z(a;,t"), with t" = nAt. Since equation (^ involves both stiff and non-stiff 
terms, it is a natural idea to employ different time-discretization strategies for 
each of them, as in 0] |2] ■ In this work we integrate ^ with a Runge-Kutta 
IMEX scheme ^^, obtaining the following semidiscrete formulation 




" v^ ' 


$2^ 


«^ 



. At 

e 



"+i = z" - At ^ fe.div/(z(^)) + — 5] hg(z^% (3) 



i=l i=l 



where the z*^*^'s are the stage values of the Runge-Kutta scheme which are given 

by 



z« = z^ 



'-' At 



At J2 a^,fcdiv/(z«) + —J2 a-fc5(^''^)> (4) 



e 

k=l fc=l 



where bi, hij and bi, aij denote the coefficients of the exphcit and iniphcit RK 
schemes, respectively. We assume that the implicit scheme is of diagonally 
implicit type. To find the z^'^'^s it is necessary in principle to solve a non linear 
system of equations which however can be easily decoupled. The system for the 
first stage z^^'^ at time t" is: 

«(i) \ / -" \ A, / \ 

u;(i) ; V ^" / ^ \ p("^'^) - ^^'^ / 

The first equation yields u^^' = u", substituting in the third equation we im- 
mediately find w^^^ and finally, substituting w^^^ in the second equation, we 
compute v^^' . In other words the system can be written in triangular form. For 
the following stage values, grouping the already computed terms in the vector 
ij(0 given by 

i-i . i-i 

SW ^ z" - At^a,,fcdiv/(z«) + _^a,,fc5(^^'^), (6) 

fe=i 

then the new stage values are given by 



e 

k=l fe=l 




= B^'^ + —au I -v^'^ + iv^e ~ D)Vw^'^ ] , (7) 

which is again a triangular system. In the numerical tests, we will apply IMEX 
schemes of order 1, 2 and 3. 

Following \2'^\ we set e = thus obtaining the so called relaxed scheme. The 
computation of the first stage reduces to 

^(i)^p(y(i)) (8) 

For the following stages the first equation is 



u(') = u''- 



At^a.^kdivv^'^l (9) 



fe=i 



In the other equations the convective terms are dominated by the source terms 
and thus w'*' and w^^' are given by 



We see that only the explicit part of the Runge-Kutta method is involved 
in the updating of the solution. Then, in the relaxed schemes we use only the 
explicit part of the tableaux. In particular we consider second and third order 
SSRK schemes ^^, namely 



IMEXl (1^* order) 



IMEX2 (2^^^ order) 



IMEX3 (3^'^ order) 
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3.1 Convergence of the semidiscrete relaxed scheme 

The aim of this section is to show the L-^ convergence of the solution of the 
semidiscrete in time relaxed scheme defined by equations ©,© and (|l()|l . We 
will extend the theorem proved in 6 , where only the case of forward Euler 
timestepping was considered. In this section, for the sake of simplicity, we set 
D = 1. 

Eliminating v from ||SIl and ^ using (|10|l . we rewrite the relaxed scheme as 









for the first stage, and 



w 









(11) 



(12) 



for subsequent stages. We recall that a Runge-Kutta scheme for the ordinary 
differential equation y' = R{y) can also be written in the form |15| 



y 



y" 



,(■0 = 



Erii«.fc(2/<'^+Aiiiti?(j/ 



(fc)- 



1 = 2,. 



.,iy, 



(13) 



where y"^"'^ = y^'^' ■ For consistency, X]fc=i '^ik = 1 for every i = 1, . . . , i/. More- 
over we assumed that aik > 0, f3ik > and that aik = implies (3ik = 0. Under 
these assumptions, each stage value y^''' can be written as a convex combination 
of forward Euler steps. This remark allows us to study the convergence of the 
Runge-Kutta scheme in terms of the convergence of the explicit forward Euler 
scheme applied to the non-linear diffusion problem. 

This latter was studied in [S] via a nonlinear semigroup argument. In the 
following we review the approach of [HI and next we extend the proof to the case 
of a (/-stages explicit Runge-Kutta scheme. 



3.1.1 The forward Euler case 

We wish to solve the evolution equation 



du 
It 



+ Lp{u) — u(-, t ~ 0) ~ uq, 



(14) 



on the domain H., where i = — A and p : R — > M is a non decreasing locally 
Lipschitz function such that p(0) — 0. Under these hypotheses, the nonlinear 
operator Au = Lpiu) with domain D{A) — {u £ L^iVt) : p{u) E D{L)} is m- 
accretive in L^(il), that is Vi^ G L^(il) and VA > there exists a unique solution 
u € DiA) such that u + \Lp{u) — ip and the application defined by (y? i-^ u is a 
contraction JJ. 

Moreover D{A) is dense in L-'^(ri), so it follows that 



SA{t)uQ = lim 

m — >oo 



-A 



Mo 



(15) 



is a contraction semigroup on L^{Vt) and 5'yi(i)uo is the generalized solution of 
H14() in the sense of Crandall-Liggett JI]. Let S{t) be the linear contraction 
semigroup generated by —L, that is u{t) — S{t)uQ is the solution of the initial 
value problem ut ~ —L{u) and u{-, t = 0) = uq. The algorithm proposed in |H] 
is 



,ri+l 



P(«") = 0, 



where t is the timestep and a^ J, 0. This can be written as 



S{<Jr) 



(16) 



,n+l _ 



= Fe{t)u" 



where FE{T)ip = ip + — [S{<7r) - l]p{^). (17) 

CTt 



Hence 



m" = (Fs(t))"uo. 



(18) 



The proof in [Sj is based on the following argument. Note that formally ^(ctt-)^ 
e~"^^ip. Let t = rn 



u{t) 



l+-±-{S{ar)-l)op 


Uq 


S'yi(Mo) when n —^ o 


Mo 





if, 



(19) 



The convergence proof requires that /i-^ < 1 where /i is the Lipschitz constant of 
p{u). We point out that ar is linked to the spatial approximation of the operator 
L and in our scheme this requirement is reflected in the stability condition of 
the fully discrete scheme (see Section 0J| 



3.1.2 Runge-Kutta schemes 



Now we are going to describe the case of a i^-stages Runge-Kutta scheme, proving 
its convergence. 

Let t > and t = t/n with n > 1; let ar ■ (0,cx)) -^ (0,oo) be a function 
such that hnv^o o^t = 0. 



^(1) — u 



Ei—l 



,(k) 



Qifc ^ 



(fc)^ 



(20) 



and proceeding as in H19f) . this becomes 



m(i) = u", 



,(») 






i('='+T§^(5K)-I)°p(«('=)) 



(21) 



We now extend ifTTI) to the Runge-Kutta scheme defined by equation l|^ . 
Define, for </> e L^{n), 

F^'Hr)<l> = El=i«.fci^('HT)0+ ^ [5(a,) -l]p(i^«(r)0), (22) 

F(r)</>==i^('')(T)</> 

and therefore 

u"it)^[FiT)ruo. (23) 

Let u(i) be the generahzed solution of (|14|l . The following theorem proves 
the convergence of the semidiscrete solution to u{t). 

Theorem 1. Assume vP £ L°°(fl), and ||m'^||oo = M ; let p be a non- decreasing 
Lipschitz continuous function on [~M, M] with Lipschitz constant /i. Assume 
that the following conditions hold 

Oiik > 0, 

P^k > 0, 

a.,k = ^ Afe = 0, 



i-l 



y aik = 1 (consistency), 



(24) 



k=l 

— < min , for r > 0, aik 7^ (stability), 

o-T (Jik 

then lim„_+oo u"(t) = u(t) in L^ . Moreover the convergence is uniform for t in 
any given bounded interval. 

The proof follows the steps of [Hj: first we show that u" verifies a maximum 
principle (Lemma ^1 and that i^ is a contraction (Lemma 121 and finally we 
apply the non linear Chernoff formula[S] . 



Lemma 1. // ^ is verified, then -M < u^ < M Vn. 

Proof. We argue by induction on n: we assume that —Ad < u^ < M and we 
show that -M < u"+^ < M. Let 

mW=fW(t)u" (25) 

Since u"+^ = u^'^^ it suffices to prove that —M < u'*^ < AI for i — 1,...,;/. 
We prove this by induction on i. When i — I, the statement is true thanks to 
the induction hypothesis on n and being F'^^^ = I. Let's assume that —M < 
y(«-i) < M holds; we are going to show that 

~M < u^''> ^ F^'\t)u" < M. (26) 

The function s i-^ aikS — ^^^^p{s) is non decreasing thanks to (|24|l and the 
hypotheses on the function p. I3y the induction hypothesis on i, we have that 
for k = I, ...,i ~ 1 



a^kM 



^p{-M) < a.,«« - ^piuC^^) < a,,M - ^p{M). (27) 



Using again the induction hypothesis on i, recalhng that p is non-decreasing, 
since 5 is a contraction in L°° 6 and p{—M) < p{u^''^) < p{M), 

p{~M) < S (p{u^''^)) < p{M) (28) 

Multiplying the last equation by ^^^ and summing it to equation \27\\ , we get 

-a,fcM<a,fewW + ^(5-I)p(7/«)<a,feM, A: = l,...,i-1 (29) 

Summing for k — 1, ...,i — 1 and using the consistency relation of H24|l: 

-M<J2 <^^ku'^''^ + —{S~ %(«('=) ) < M (30) 

fc=i '^^ 

In particular this is valid when i = v, proving that —M < m("+-'^) < M. D 

Now we can replace p by p, where p = p in — M < x < M, p — p{M) for 
X > M and p — p{—M) for x < —M : the algorithm is the same and in what 
follows we can assume that p is Lipschitz continuous with constant ^ on all M. 

Lemma 2. If the hypotheses of Theorem^ hold, then F{t) is a contraction on 
L^{n), I.e. 

||F(T)0-F(r)V||i<||0-Vl|i V7/;,0eii (31) 
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Proof. We start showing that the result holds for a single forward Eiiler step. 
Recalling the definition of Fe from (|17|l 

\\Fe{t)^^Fe{t)^i < ^ \\S{ar)[pW - pimii + ||(0- V') - ^[p{^)~pm 

< ^ \\p{<i>) - pinii + \\{<t'- ^pw) - (^ - ^pw) II ^ 

= ll0-^lli 

(32) 
where we used the contractivity of S. The last equality relies on the fact that p 
and the function x i— > x — —p{x) are non-decreasing, which in turn is guaranteed 
by the stability condition, that in this case reduces to nr/ar < 1 [HI- 
In the general case we have: 

Pe (iff) F^'Hr)cly~FE (^) F^'^Ht)^ 



<U-i'\\i 

(33) 
In the second inequality we used the contractivity of Fe and the stability condi- 
tion, while in the third one we apply an induction argument on the contractivity 
of F^''' , the positivity constraint on aik and /3ik, as well as the consistency con- 
dition J2k ^ik = 1- Setting i = v yields the result. D 

Proof of Theorem^ Let f/'r and V' be respectively 



■01- 



/+-(/-F(r)) 



and 



^ = {I + \Ay^ (j). 



(34) 



The function ijj exists since the operator A is m-accretive, whereas the existence 
of the function ipT is guaranteed by the following fixed-point argument. Let 



G{v) 



1 



V 



l+lf !]+! 



F{T)y 



where (j) £ L^, y E D{A) and rj > 0. We have. 



\\G{y)-G{x)\\ = 



V 



\F{T)y - FiT)x\\ < 



V 



■\\y-^\\ 



ri + I 1] + I 

since i^ is a contraction, as proved in Lemma El Thus G is also a contraction 
and therefore it possesses a unique fixed point which coincides with ■(/;,-. 
We want to show that 

■0^ — > ^ in L 

as r -^ for each fixed A > 0. Let 

(br^^+-il-F{T)),p. 
T 

We want to estimate ^t- — i/; in terms of 0r — 0- 

0r - = (1 + -)(V^ - A) - -{F{T)i, - F{T)i,,) 
T T 
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Therefore 



(l + -)(V'-^r) 



^)^-{F{t)^-F{t)^,) 



and taking norms and using the fact that F is contraction we have 



(l + -)ll^-V^r||-||</'r- 



<||(i + -)(^-V.) 

T 



m<-u~4'r\\ 

T 



In particular 



(l + -)IIV'-V^r||-||<^r-<^||<-||V'-^r|l 



and therefore 1 1 ^ — V-'t 1 1 < 
Now we estimate ||0 — 
Note that 



in the simple case of a forward Euler scheme. 



A, 



^,^XAi;^-{l-FiT))^ 



and thus \\(f> — 0r|| measures a sort of consistency error. For a single forward 
Euler step, F = Fe where Fe is defined in ifTTj) . Thus 



= A 



A,jj-—{l-S{ar))p{^) 







(35) 



as r ^ since 



I-Sja^) 



'P{ip) -^ Lp{il)) — A'ip. 



The more general case of a z^-stages Runge-Kutta scheme can be carried 
out by induction following the procedure already applied in the proofs of the 
previous lemmas. 

We now use Theorem 3.2 of jHj which, specialized to our case, can be written 
as follows. Assume that F{t) : L^ ^ L^ for r > is a family of contractions. 
Assume further that an m-accretive operator A is given and let S(t) be the semi- 
group generated by A. Assume further that the family F{t) and the operator 
A are linked by the following formula 



for each 



e Li. Then 



A , , , , 

I + -{I -Fir)) 



(f, ^ ^Jj ^ (I + XAY 



(36) 



jim^ F{-] (j) = S{t)(t) 



'&L\ 



D 



4 Fully discrete relaxed scheme 

In order to complete the description of the scheme, we need to specify the space 
discretization. We will use discretizations based on finite differences, in order 
to avoid cell coupling due to the source terms. 
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Note that the IMEX technique reduces the integration to a cascade of re- 
laxation and transport steps. The former are the impHcit parts of lO and (O, 
while the transport steps appear in the evaluation of the explicit terms B^^' in 
(|HJl. Since jSJ) and Q involve only local operations, the main task of the space 
discretization is the evaluation of div(/), where we will exploit the linearity of 
/ in its arguments. 

4.1 One dimensional scheme 

Let us introduce a uniform grid on [a, b] G M., Xj = a — ^ + jh for j — 1, . . . ,n, 
where h = {b — a)/n is the grid spacing and n the number of cells. The fully 
discrete scheme may be written as 

i=l i=l 

where F , ,„ are the numerical fluxes, which are the only item that we still need 

J + l/2 ' ■' 

to specify. For convergence it is necessary to write the scheme in conservation 
form. Thus, following |33], we introduce the function F such that 



f[z{x,t))^-j F{s,t)As ^^{z{xj,t))^-i^F{xj+,^2,t)-F{xj_,/^,t)j 



The numerical flux function Fj^i/2 approximates F{xj^i^2)- 

In order to compute the numerical fluxes, for each stage value, we reconstruct 
boundary extrapolated data z L/2 with a non-oscillatory interpolation method 

from the point values z^^' of the variables at the center of the cells. Next we 
apply a monotone numerical flux to these boundary extrapolated data. 

To minimize numerical viscosity we choose the Godunov flux, which in the 
present case of a linear system of equations reduces to the upwind flux. In 
order to select the upwind direction we write the system in characteristic form. 
The characteristic variables relative to the eigenvalues (p,—(p,Q (in one space 
dimension if reduces to a scalar parameter) are respectively 

V + ipW WW — V 

U = ^ V ^ W ^u-w. (2) 

Note that u = U + V + W. Therefore the numerical flux in characteristic 
variables is -P,+i/2 = (<^[r^^/2' ^'^^J>l/2'0)■ 
The accuracy of the scheme depends on the accuracy of the reconstruction 
of the boundary extrapolated data. For a first order scheme we use a piecewise 
constant reconstruction such that U^,^,^ — Uj and V .^j^ — Vj+i- For higher 
order schemes, we use ENO or WENO reconstructions of appropriate accuracy 

(123). 

For e — > we obtain the relaxed scheme. Recall from equation (|10|l that the 
relaxation steps reduce to 

wf=p{uf), vf^~DVwf, (3) 

13 



where V is a suitable approximation of the one-dimensional gradient operator. 
Thus the transport steps need to be applied only to u'-'-' 



-f - < 



- A^a.. ^ 0^:^, - U^-^ - . {V^-, - l^^i), (4) 



fe=i 
Finally, taking the last stage value and going back to conservative variables, 

,n+l 



"j 



+ '^K + l/2 - ^J + 1/2 - (^J-1/2 - %-l/2)]J 

We wish to emphasize that the scheme reduces to the time advancement 
of the single variable u. Although the scheme is based on a system of three 
equations, the construction is used only to select the correct upwinding for 
the fluxes of the relaxed scheme and the computational cost of each time step 
remains moderate. 

4.2 Non linear stability for the first order scheme 

The relaxed scheme in the first order case reduces to: 

(6) 
We wish to compute the restrictions on A and (p so that the scheme is total 
variation non-increasing. We select the centered finite difference formula to 
approximate the partial derivatives of p(u); we drop the index n and write pj 
for p{u'[j). Define Aj_|_i/2 = ^^.^^_^^. and observe that these quantities are always 
nonnegative since p is nondecreasing. We obtain 

+ (i-^(27r + '^)^j-i/2)l%-%-il 

(7) 
provided that 

1 

'2h 

Assuming that the data have compact support, we can rescale all sums and 
finally get TV{u"~^^) < TV(m"). Taking into account the Lipschitz condition 
on p, the scheme is total variation stable provided that © is satisfied, i.e. that 

2/i2 1 (2 - (5) o 

At < -— ~ -^ ^-h^ (9) 

~ fi 1 + 2hip fi ^ ' 

where 5 vanishes as h does. We point out that the stability condition is of 
parabolic type. Finally, we observe that using one-sided approximations for the 
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(1_A(— +<^)A,_i/2)>0 Vj (8) 



partial derivatives of p in the scheme ©, one gets a stabihty condition involving 
the relation ip > 1/h. This would reintroduce in the scheme the constraint due 
to the stiffness in the convective term that prompted the introduction of ip in 

4.3 Linear stability 

We study the linear stability of the schemes based on equations , Ql and lO) 
in the case when p{u) — u, by von Neumann analysis. We substitute the discrete 
Fourier modes u" = p^Q^Uf^/^) into the scheme, where k is the wave number and 
N the number of cells. Let ^ — k/N, we compute the amplification factor Z{£^) 
such that u"^^ — Z(^)w". We can consider ^ as a continuous variable, since 
the amplification factors for various choices of N all lie on the curves obtained 
considering the variable ^ € [0, 27r]. 

First we consider the same scheme studied in the previous section, for com- 
parison purposes. Using piecewise constant reconstructions in space and forward 
Euler time integration, the amplification factor is Z{£^) — 1 + M(^), where 

M(0 - ^ (cos(0 - 1) (cos(0 + 1 + hp) . 

M{^) takes maximum value and attains its minimum at the point ^* such that 
cos(^*) — ~iph/2. Stability requires that M(^*) > —2, i.e. 

and recalling that A = At/ ft,, 

2/)2 
At < 2 - 2 (1 - ph) /i2 (10) 



(l + f)' 



This gives a CFL condition of the form At < 2(1 — S)h'^ where 6 = 0{h(p) (see 
figure ^). These results are in very good agreement with those of the nonlinear 
analysis performed in the previous section. 

Now we consider higher order spatial reconstructions coupled with forward 
Euler timestepping. M takes the form 

A^(e,7) = ^[/i(cos(0)+7/2(cos(C))] 

where 7 = hip. Since 7 is small, we compute the critical points ^* of Af(^,0). 
For stability we thus require that —2 < M{^* , 7) < 0. 

We consider a piecewise linear and a WENO reconstruction. The first one is 
computed along characteristic variables using the upwind slope, while the gra- 
dient oi p{u) is computed with centered differences. The WENO reconstruction 
is fifth order accurate and is obtained by setting to 1 the smoothness indicators 



15 




Figure 1: Amplification factor for upwind spatial reconstruction coupled with 
forward Euler (left) and for upwind second order spatial reconstruction coupled 
with second order time integration (right). 

and the gradient of p{u) is computed with the fourth order centered difference 
formula. 

For the piecewise linear reconstruction, we have that 



A^(0 = -j; [(cos2(e) - l)(cos(0 - 2) + h^icosiO - If 

and therefore 

2/i2 



At < 



20+14^/7 I 8+2^7. 



0.94(1 - 1.44(/?/i)/i^ 



27 +'-^f^^h 

For the WENO reconstruction M{^,j) can be easily computed and we get 

At < 0.79(1 - 0.13tph)h'^ 

Now we wish to extend our results to the case of higher order Runge-Kutta 
schemes. Since both the equation and the scheme are linear, the amplification 
factors for the Runge-Kutta schemes of order 2 and 3 used here are respectively 

^(2)(e) = i + M(0 + ^ 

where M(^) is the function appearing in the amplification factor relevant to the 
chosen spatial reconstruction. We have that 

Z('2)(e)-M'(e)(i + M(e)) 

and therefore the critical points are the points £,* such that M'(^*) = 0. 
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/ 


^\ '7 V 
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Figure 2: Amplification factors Z for WENO reconstructions of order 5 coupled 
with first order (left) and third order (right) time integration. 



In the Runge-Kutta 2 case the stability constraint ||.Z'(2)(^*)|| < 1 reduces to 
the CFL condition for the forward Euler scheme. For Runge-Kutta 3, |j^(3)('f*)|| < 
1, provided that 

Af(^*) > S~ -2.51 

Notice that this is less restrictive than the Euler and RK2 schemes for which 
the stability requirement is M(^*) > —2. 

For the Runge-Kutta 3 scheme with linearized WENO of order 5, we have 



At < 



2.51 + 0.33v3/i 



(1 - .Vyihiph)h^ 



Table n summarizes the stability results obtained in this section listing the 
values of the constant C that appear in the stability restriction At < C(l — 
Ci(ph)h^. Figures n and 121 contain the amplification factors Z{£^) for tp = 1 
and h = 10~^ for various choices of spatial reconstructions and time integration 
schemes. Each of them contains the curve corresponding to the value of C 
reported in Tableland two other close- by values. 





RKl 


RK2 


RK3 


P-wisc constant 


2 


2 


2.51 


P-wisc linear 


0.94 


0.94 




WEN05 


0.79 


0.79 


1 



Table 1: 



4.4 Boundary conditions 

Different boundary conditions can be implemented. Here we describe how to 
implement Neumann boundary conditions, considering for simplicity the one- 
dimensional case. 
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We first add g ghost points on each side of the computational domain [a, b], 
where g depends on the order of the spatial reconstruction. We find a polyno- 
mial q{x) of degree d passing through the points {xi, Ui) for i = 1, . . . , d and 
having prescribed derivative at the boundary point a;i/2 = a- (The degree d is 
determined by the accuracy of the scheme that one wants to obtain and should 
match the degree of the reconstruction procedures used to obtain U^ and V^ .) 
This polynomial is then used to set the values U-i = q{x-i) of the ghost points 
for i = 0, 1, (7— 1. One operates similarly at the right edge of the computational 
domain. 

We also used periodic boundary conditions, which can be implemented with 
obvious choice of the values Ui at the ghost points. 

4.5 Multi-dimensional scheme 

An appropriate numerical approximation of Q in K'^ that generalizes the scheme 
described in Section Wl\ can be obtained by additive dimensional splitting. We 
consider the relaxed scheme, i.e. e = and for the sake of simplicity, let us 
focus on the square domain [a, b] x [a, b] C M'^. Here we shall describe the gen- 
eralization of the scheme defined by equations Q, @, Q and jSJ to the case 
of two space dimensions. 

Without loss of generality, we consider a uniform grid in [a, 5] x [a, 6] C M^ 
such that Xij = {xi,yj) = {a — h/2,a-'h/2)+i{h,0)+j{0, h) ior i,j ^ 1,2, ... ,n 
and h = {b — a)/n. 

In the present case, u and w are one-dimensional variables, while v = 
{v{i),V(2)) is now a field in R^. First we observe that the relaxation steps lj3J) 
are easily generalized for d > 1. For the transport steps, one has to evolve in 
time the system 





( u \ 




d 


V(l) 


d 

-\ — 


dt 


V(2) 
\ w J 


ox 



10 

(/^^ 



10 





( u \ 






«(1) 


d 




V(2) 


'^dy 




\ ^ J 





10 



^2 

10 





( u \ 




W(l) 




V(2) 




\ w J 



The semidiscretization in space of the above equation can be written as 



(11) 



= 



dz. 



"i-,3 



1 



1 



dt 



{Fi+l/2.j - F^_i/2^j) - - (G'ij + 1/2 - G.,j_i/2) 



where F and G are the numerical fluxes in the x and y direction respectively 
and can be written as 



F,_ 



-l/2j 



F(z 



i+l/2j"' i+l/2,jl 



G 



■'j+1/2 



G{z 



+ z' ) 

ij + l/2' i,j + l/2' 



The fluxes in the two directions are computed separately. We illustrate the 
computation of the flux F along the x direction. We note that only the fleld V(i) 
appears in the differential operator along this direction. The third component 
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of the flux is zero and thus we have three independent characteristic variables, 
namely 

LpW + Vi (pW-Vi 

which correspond respectively to the eigenvalues 1^9, —1^9,0. At this point the 
numerical fluxes can be easily evaluated by upwinding. We proceed similarly for 
the numerical flux G that depends on the characteristic variables f/(2), ^{2), W. 



Xi 






Denote by ^7^^2/27 *^^ reconstructions of C7(i)(-,2/j) at the point ( 
h/2,yj). This involves a reconstruction of the restriction of t/(i) to the line 
y = Hi and can be obtained with any of the one-dimensional techniques men- 
tioned in Secti0n f4.ll Similarly, denote U^ +1/2 ^^'^ reconstructions of U(^2){xi, ■) 
at the point {xi,yj + h/2). Now, formulas Q and lO become respectively 

'P\^^,] + l/2 '^t,j-l/2) '^l^^j + 1/2 ^J-1/2J 

(12) 
and 

"ij - "*.i " -^ 2^1=1 '■P°l [\}'i+l/2,j - ^+l/2j J - y^i-1/2,3 - ^i-l/2,3) 

(rjH)- _y{l)+ \_(uW- _yW+ \ 
\'^zj + l/2 ^i,J + l/2j V JJ-1/2 ^ij-1/2) 

(13) 
The generalization to d > 2 and rectangular domains is now trivial. We 
stress once again that no two-dimensional reconstruction is used, but only d one- 
dimensional reconstructions are needed. Finally, boundary conditions can be im- 
plemented direction-wise with the same techniques used in the one-dimensional 
case. 

5 Numerical results 

We performed several numerical tests of our relaxed schemes. First we tested 
convergence for a linear diffusion equation with periodic and Neumann boundary 
conditions for initial data giving rise to smooth solutions. Next, numerical tests 
were also performed on the porous media equation ut = {u™)xx, m = 2, 3, both 
in one and two dimensions. 

5.1 Linear diffusion 

For the first test we considered the linear problem 

du , . d^u , , .„ , 

-^{x,t) ^ —u{x,t) xe[0,l] 

u{x, 0) = uq{x) X G [0, 1] 

First we used periodic boundary conditions with U{){x) = cos{2'itx), so that 
u{x, i) = cos.{2i:x)e^'^'^ * is an exact solution. Then we used Neumann boundary 
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N=40 


N=80 


N=160 


N=320 


N=640 


EN02, RKl 


2.012e-03 


5.6378e-04 


1.0736e-04 


1.5539e-05 


2.5065e-06 


EN03, RK2 


1.9066e-06 


2.3057O-07 


5.6115C-08 


8.6904e-09 


1.1905e-09 


EN04, RK2 


7.7517C-06 


5.7082C-07 


3.3507C-08 


1.4978e-09 


7.0725e-ll 


EN05, RK3 


1.3864e-08 


6.0259C-10 


2.2121C-11 


7.4454e-13 


2.3803e-14 


EN06, RK3 


1.5538e-08 


8.5661C-10 


1.4460-11 


1.7111e-13 


1.5311e-15 


WEN03, RK2 


1.9799C-03 


5.1278c0-4 


1.4332e-04 


2.1488e-05 


7.512e-08 


WEN05, RK3 


1.5892e-07 


4.8069e-09 


1.59e-10 


5.2337e-12 


1.6758e-13 





N=40 


N=80 


N=160 


N=320 


N=640 


EN02, RKl 


1.3973 


1.8354 


2.3926 


2.7886 


2.6322 


EN03, RK2 


5.9501 


3.0477 


2.0388 


2.6909 


2.8678 


EN04, RK2 


3.8987 


3.7634 


4.0905 


4.4836 


4.4045 


EN05, RK3 


6.8124 


4.524 


4.7677 


4.8929 


4.9671 


EN06, RK3 


5.9907 


4.181 


5.8885 


6.401 


6.8043 


WEN03, RK2 


0.56648 


1.949 


1.8391 


2.7376 


8.1601 


WEN05, RK3 


2.9595 


5.0471 


4.918 


4.925 


4.9649 



Table 1: L^ norms of the error and convergence rates for the hnear diffusion 
equation with periodic boundary conditions, with smooth initial data 

conditions Ux{0) = Ux{l) ~ 1 with initial data uo{x) ^ x + cos{2t:x), so that 
u{x,t) = X + cos{2TTx)e~'^^ * is an exact solution. 

We tested the numerical schemes defined by equations (0), ®j Q) and 
(j21) with various degrees of accuracy for the spatial reconstructions and time- 
stepping operators. We used ENO spatial reconstructions of degrees from 2 to 
6 and WENO reconstructions of degrees 3 and 5. The time-stepping procedures 
chosen are IMEX Runge-Kutta schemes of SectionOof accuracy chosen to match 
the accuracy of the spatial reconstruction. Since stability forces the parabolic 
restriction At < Ch? , an IMEX scheme of order m was coupled with a spatial 
ENO/WENO reconstruction of accuracy p such thatp < 2m, obtaining a scheme 
of order p. 

We computed the numerical solution of the diffusion equation with final time 
t = 0.05 with N = 40, 80, 160, 320, 640 grid points and computed the L^ norm of 
the difference between the numerical and the exact solution. The results are in 
Tables n for the periodic boundary conditions and 12 for the Neumann boundary 
conditions. One can see that the expected convergence rates are reached, even 
if the combination of the WENO reconstruction of accuracy 3 and the IMEX 
scheme of second order reach the predicted error reduction only on very fine 
grids. 

5.2 Porous media equation 

On the porous media equation (^ with p{u) = u™ we performed a test proposed 
in ^j. We took m = 2, 3 and initial data of class C^ as follows: 



m(x,0) = 



cos^(7ra;/2) 




|x|<l 
\x\ >1 



(1) 
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N=40 


N=80 


N=160 


N=320 


N=640 


EN02, RKl 


2.1965e-03 


5.7152e-04 


1.4301e-04 


2.32e-05 


4.743e-06 


EN03, RK2 


2.0621C-06 


2.2641C-07 


6.7935e-08 


8.8255O-09 


1.2339C-09 


EN04, RK2 


8.1764C-06 


5.4431C-07 


3.6974e-08 


1.3686C-09 


8.335e-ll 


EN05, RK3 


1.5484C-07 


4.4163C-09 


1.2405e-10 


3.7803O-12 


1.1669C-13 


WEN03, RK2 


1.9092C-03 


4.4225e-04 


1.2914e-04 


4.5037O-06 


7.4526e-08 


WEN05, RK3 


2.5048e-07 


4.9279e-09 


1.4776e-10 


4.7482e-12 


1.4948e-13 





N=40 


N=80 


N=160 


N=320 


N=640 


EN02, RKl 


1.4361 


1.9424 


1.9987 


2.624 


2.2902 


EN03, RK2 


6.1004 


3.1871 


1.7367 


2.9444 


2.8385 


EN04, RK2 


3.9763 


3.909 


3.8798 


4.7558 


4.0373 


EN05, RK3 


5.6626 


5.1317 


5.1539 


5.0362 


5.0178 


WEN03, RK2 


1.2624 


2.11 


1.7759 


4.8417 


5.9172 


WEN05, RK3 


4.9122 


5.6676 


5.0597 


4.9597 


4.9893 



Table 2: L^ norms of the error and convergence rates for the hnear diffusion 
equation with Neumann boundary conditions, with smooth initial data. 

The computational domain is {|a;| < 3} C R and the boundary conditions are 
periodic; the CFL constant is taken as C = 0.25. 

Since the initial data has compact support and is Lipschitz continuous, the 
solution will be of compact support for every t > 0, but will develop a disconti- 
nuity in Ux at some finite time r > (see PD- 
As was shown in j3), the solution with the initial condition we chose has 
a front that does not move for t < 0.034. We therefore chose a final time 
of the simulation tan = 0.03 to prevent the formation of the singularity of Ux 
from affecting the order of convergence. We used as reference solution the one 
obtained numerically with TV = 4860 grid points and computed the L^ norms of 
the errors of the solutions with TV = 60, 180, 540, 1620 grid points. The results 
are presented in Tabic 13 

First of all one verifies that the degree of regularity of the solution poses 
a limit on the order of convergence of the schemes: therefore the schemes we 
tested perform at best as third order schemes, as confirmed by the data in Table 
El Still, high order schemes yield smaller error on a given grid. This can be 
of practical importance in problems where one does not have the freedom of 
choosing the number of grid points, as in digital image analysis, where non- 
linear degenerate diffusion equations are sometimes used as filters for contour 
enhancement (see jH])- 

In Figure n we show the numerical solution for the porous media equation 



with p{u) = u and p{u) 



with the initial data ^ and t £ [0, 2]. It can 



be appreciated that a front (i.e. a discontinuity of ^) develops at a finite time 
and then it travels at finite speed. 

We present a numerical simulation for the two-dimensional porous media 
equation QJ with p{u) = v? . We chose an initial data uo{x,y) given by two 
bumps with periodic boundary conditions on [—10,10] x [—10,10]. The large 
domain ensures that the compact support of the solution is still contained in 
the computational domain at the final time of the calculation. The numerical 
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N=60 


N=180 


N=540 


N=1620 


EN02, RKl 


2.6365e-04 


1.9898e-05 


2.049e-06 


2.076e-07 


EN03, RK2 


1.9605e-05 


6.0423C-07 


2.4141C-08 


8.9729e-10 


EN04, RK2 


1.2127e-05 


2.967C-07 


9.9925C-09 


3.5781e-10 


EN05, RK3 


4.694e-06 


1.719O-07 


6.3248O-09 


2.4447e-10 


EN06, RK3 


4.1099e-06 


1.4711e-07 


5.3992e-09 


2.0849e-10 


WEN03, RK2 


1.5871e-04 


1.0448e-05 


4.3463e-07 


8.8767e-09 


WEN05, RK3 


7.5662e-06 


4.6049e-07 


7.4746e-09 


2.7985C-10 





N=60 


N=180 


N=540 


N=1620 


EN02, RKl 


2.8243 


2.352 


2.0692 


2.084 


EN03, RK2 


5.1899 


3.1672 


2.931 


2.9968 


EN04, RK2 


5.6271 


3.3774 


3.0865 


3.0307 


EN05, RK3 


6.491 


3.0103 


3.006 


2.9611 


EN06, RK3 


6.612 


3.0311 


3.0083 


2.962 


WEN03, RK2 


3.2863 


2.4765 


2.8942 


3.5418 


WEN05, RK3 


6.0565 


2.5479 


3.7509 


2.9902 



Table 3: L^ norms of the error and convergence rates for the porous media 
equation periodic boundary conditions, with initial data of class C^. 





0.2 



Figure 1: Snapshots of the numerical solutions for the porous media equation 
with p{u) = V? (left) and p{u) = v? (right). Initial data are chosen according 
to (0 and the numerical solutions are represented at times i = 0, 0.2, . . . , 2.0. 
The solutions are obtained with the spatial WENO reconstruction of order 5 
and the RK3 time integrator. 
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Figure 2: The numerical solution of the porous media equation on a square 
regular grid with compactly supported initial data. From top left to bottom 
right, we show the numerical solution at times t := 0, 0.5, 1.0, 4.0. 
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approximation at different time levels is shown in Figure |21 

We can note that the symmetries of the initial data are preserved and the so- 
lution seems to be unaffected by the dimensional splitting of the two-dimensional 
scheme. 



6 Conclusions 

We have proposed and analyzed relaxed schemes for nonlinear degenerate parabolic 
equations. 

By using suitable discretization in space and time, namely ENO/WENO 
non-oscillatory reconstructions for numerical fluxes and IMEX Runge-Kutta 
schemes for time integration, we have obtained a class of high order schemes. We 
have developed a theoretical convergence analysis for the semidiscrete scheme; 
furthermore we studied stability for the fully discrete schemes. Our computa- 
tional results suggest that our schemes converge with the predicted rate. 

Finally, we point out that these schemes can be easily implemented on par- 
allel computers. Some preliminary results and details are reported in [S]. In 
particular the schemes involve only linear matrix-vector operations and the ex- 
ecution time scales linearly when increasing the number of processors. 

Our numerical approach can be easily extended also to more general prob- 
lems, as nonlinear convection-diffusion equations or nonlinear parabolic systems. 
Some of these applications will appear in a forthcoming paper. 
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